{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "35302d04",
   "metadata": {},
   "source": [
    "# test 251119"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc0e141b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from icecream import ic\n",
    "from pathlib import Path\n",
    "\n",
    "_CWD = Path().cwd()\n",
    "_PRJROOT = _CWD.parent.parent\n",
    "_G4OUT = _PRJROOT / \"G4\" / \"build\" / \"Release\"\n",
    "\n",
    "data_path = _G4OUT / \"LDose_Z2_275.00MeV.npz\"\n",
    "linear_data = np.load(data_path)\n",
    "## Z bin; Phi bin; R Bin; 右边的指标变化最快\n",
    "ix0, ix1, ix2, val, val2, entry = linear_data[\"idx0\"], linear_data[\"idx1\"], linear_data[\"idx2\"], linear_data[\"value\"], linear_data[\"val2\"], linear_data[\"entry\"]\n",
    "\n",
    "## test\n",
    "if False:\n",
    "    ic(ix0[:12])\n",
    "    ic(ix1[:12])\n",
    "    ic(ix2[:12])\n",
    "    ic(val[:12])\n",
    "\n",
    "##xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\n",
    "Zbin: int = 600\n",
    "Zmax: float = 300\n",
    "Zmin: float = -300\n",
    "## 在 radial 方向, 分成 200 layers; 角度方向不分sector, 只有一整圈\n",
    "Rbin: int = 200\n",
    "Rmax: float = 200\n",
    "Rmin: float = 0\n",
    "##\n",
    "Rmax_fit: float = Rmax\n",
    "Rmin_fit: float = Rmin\n",
    "###=========================================\n",
    "ZDiff: float = Zmax - Zmin\n",
    "Zbwth: float = ZDiff / Zbin  # z bin width\n",
    "Zaxis = np.arange(Zmin + Zbwth / 2.0, Zmax, Zbwth)\n",
    "###=========================================\n",
    "RDiff: float = Rmax - Rmin\n",
    "Rbwth: float = RDiff / Rbin\n",
    "Raxis = np.arange(Rmin + Rbwth / 2.0, Rmax, Rbwth)\n",
    "\n",
    "###=========================================\n",
    "material = \"WATER\"\n",
    "prtcls_name = \"he\"\n",
    "# enedir = '125MeV'\n",
    "dir_tag = \"filter\"\n",
    "\n",
    "ic(val.shape)\n",
    "ic(val.flags[\"C_CONTIGUOUS\"])\n",
    "## G4 dump 结果排布次序为 C order, last fastest; # iZ, iPHI, iR\n",
    "data_2D = val.reshape((Zbin, Rbin), order=\"C\")\n",
    "data_intgrd = data_2D.sum(axis=1)  ## dicrete integrad\n",
    "## inspect data\n",
    "ic(data_2D.shape)\n",
    "ic(data_2D.flags[\"C_CONTIGUOUS\"])\n",
    "ic(data_intgrd.shape)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(data_intgrd)\n",
    "ax.set_xlim(0, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccaf10ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "if False:\n",
    "    import importlib\n",
    "    import pyauto.fit.fitmodel\n",
    "    import pyauto.fit.multigauss\n",
    "    import pyauto.fit.fitutils\n",
    "\n",
    "    importlib.reload(pyauto.fit.multigauss)\n",
    "    importlib.reload(pyauto.fit.fitmodel)\n",
    "    importlib.reload(pyauto.fit.fitutils)\n",
    "\n",
    "from pyauto.fit.multigauss import fit_multi_gauss, multi_gauss_model\n",
    "from pyauto.fit.fitmodel import (\n",
    "    FitConfig,\n",
    "    RunMGFitIn,\n",
    "    MGaussIn,\n",
    "    myDftVals as mydv,\n",
    ")\n",
    "from pyauto.fit.fitutils import (\n",
    "    normalize_by_sum,\n",
    "    NmlzIn,\n",
    "    MeshEnum,\n",
    ")\n",
    "\n",
    "# 剂量最大位置\n",
    "idz0 = np.argmax(data_intgrd)\n",
    "peak_idd = data_intgrd[idz0]\n",
    "ic(idz0)\n",
    "ic(peak_idd)\n",
    "\n",
    "\n",
    "z_depth_index = int(0.85 * idz0)  # 选一个peak之前的位置\n",
    "print(\"Zaxis[{}]: {}\".format(z_depth_index, Zaxis[z_depth_index]))\n",
    "\n",
    "\n",
    "# fn_normalize = normalize_by_sum\n",
    "# 归一化数据\n",
    "profile = data_2D[z_depth_index, :]\n",
    "profile = normalize_by_sum(NmlzIn(mtype=MeshEnum.Cyliner, vars_tgt=profile, bin_width=Rbwth))\n",
    "# ic(profile)\n",
    "##=====================================================\n",
    "## 限制拟合数据的横坐标范围\n",
    "mask = (Raxis >= Rmin_fit) & (Raxis <= Rmax_fit)\n",
    "Raxis_fit = Raxis[mask]\n",
    "profile = profile[mask]\n",
    "\n",
    "\n",
    "##=====================================================\n",
    "# N-高斯拟合设置\n",
    "n_gauss = 5\n",
    "# fit results 存储所有深度的拟合结果\n",
    "fit_rslts = np.zeros((Zbin, 2 * n_gauss + 1))\n",
    "output_dir = Path(\"fitresults\")\n",
    "output_dir.mkdir(exist_ok=True, parents=True)\n",
    "## 先拟合布拉格峰位附近\n",
    "crds, vals_tgt = Raxis, profile\n",
    "# fitData = FitData(crds, vals)\n",
    "# fitConfig = FitConfig(n_gauss)\n",
    "fitCfg = FitConfig(\n",
    "    ngauss=5,\n",
    "    with_background=False,\n",
    "    is_post_bragg=False,\n",
    "    lam_n=mydv.penalty_norm,\n",
    "    info_verbose=True,\n",
    "    info_tag='11',\n",
    ")\n",
    "fitRet = fit_multi_gauss(\n",
    "    RunMGFitIn(\n",
    "        crds=crds,\n",
    "        vals_tgt=vals_tgt,\n",
    "        fitcfg=fitCfg,\n",
    "        sigma_prev=None,\n",
    "        gpars_init=None,\n",
    "    )\n",
    ")\n",
    "\n",
    "gpars = fitRet.gpars2\n",
    "# ic(gpars)\n",
    "print(\n",
    "    \"==============================================\\n\",\n",
    "    f\"normFactor: {gpars.norm_factor}\\n\",\n",
    "    f\"第一高斯权重: {gpars.w0}\\n\",\n",
    "    f\"其余权重   : {gpars.weights}\\n\",\n",
    "    f\"其余权重和 : {gpars.sum_weights2n}\\n\",\n",
    "    f\"sigmas   : {gpars.sigmas}\\n\",\n",
    ")\n",
    "vals_p = multi_gauss_model(MGaussIn(crds=crds, gpars=gpars))\n",
    "\n",
    "\n",
    "y_min = np.min(vals_tgt)\n",
    "y_max = np.max(vals_tgt)\n",
    "y_lower = max(y_min * 0.5, 1e-12)\n",
    "y_upper = y_max * 2\n",
    "plt.figure(dpi=150)\n",
    "plt.semilogy(crds, vals_tgt, \"o\", label=\"data\", markersize=4, markerfacecolor=\"none\", alpha=0.8)\n",
    "plt.semilogy(crds, vals_p, \"-\", label=\"fit\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"x [mm]\")\n",
    "plt.ylabel(\"Dose [a.u.]\")\n",
    "plt.title(f\"Depth idx={z_depth_index} (z={Zaxis[z_depth_index]:.3f} cm)\")\n",
    "plt.ylim(y_lower, y_upper)\n",
    "plt.show()\n",
    "plt.savefig(f\"{output_dir}/depth_{z_depth_index:03d}.png\")\n",
    "# plt.savefig(f\"depth_{z_depth_index:03d}.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c0e1d3d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyauto",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
